#### "Inclusionary regimes, party institutionalization and redistribution under authoritarianism" ####
# authors: "Lars Pelke"
# date: 2020-06-16
# written under "R version 3.6.1 (2019-07-05)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

#libraries
library(tidyverse)
library(ggpubr)
library(texreg)
library(dotwhisker)
library(broom.mixed)
library(countrycode)
library(readstata13)
library(imputeTS)


library(lme4)
library(sjPlot)
library(sjmisc)
library(lfe)
library(ggeffects)
library(stargazer)


# set working directory
# please use the working directory, where you stored the zip-file. 

##########################################################################
##########################################################################

#### Outline ####

## 1. Small Sample based on Solt 2019 original data
## 2. REWB with democratic and autocratic countries 
## 3. Income Taxation Data 
## 4. Social Spending Data 
## 5. Social Protection Data
## 6. Public Goods Provision V-Dem Data
## 7. OLS Models relative redistribution

##########################################################################
##########################################################################

#### APPENDIX B ####
#### 1. Small Sample based on Solt 2019 original data ####

## Load Data ##

## V-DEM
vdem <- read_csv("data/vdem_9/V-Dem-CY-Full+Others-v9.csv") 

vdem$cown <- countrycode(vdem$country_name, "country.name", "cown", warn = TRUE)
vdem$cown[vdem$country_name == "Hong Kong"] <- "715"
vdem$cown[vdem$country_name == "Republic of Vietnam"] <- "817"
vdem$cown <- as.integer(vdem$cown)

## Solt Inequality and Redistribution Data 

swiid <- load("data/Solt 2019 SWIID 8.2/swiid8_2.rda")
swiid8.0 <- read_csv("data/Solt 2019 SWIID 8.0//swiid.csv") 

swiid_summary$country_year <- paste(swiid_summary$country, swiid_summary$year, sep = "")
swiid8.0$country_year <- paste(swiid8.0$country, swiid8.0$year, sep = "")

setdiff(swiid_summary$country_year, swiid8.0$country_year)

swiid_summary$cown <- countrycode(swiid_summary$country, "country.name", "cown", warn = TRUE)
swiid_summary$cown[swiid_summary$country == "Hong Kong"] <- "715"
swiid_summary$cown[swiid_summary$country == "Puerto Rico"] <- "1014"
swiid_summary$cown[swiid_summary$country == "Anguilla"] <- "1022"
swiid_summary$cown[swiid_summary$country == "Micronesia"] <- "987"
swiid_summary$cown[swiid_summary$country == "Palestinian Territories"] <- "1020"
swiid_summary$cown[swiid_summary$country == "Serbia"] <- "345"
swiid_summary$cown[swiid_summary$country == "Turks and Caicos Island"] <- "1018"
swiid_summary$cown <- as.integer(swiid_summary$cown)


summary(swiid_summary$rel_red)
summary(swiid_summary$abs_red)

## Correlates of War National Material Capabilities Dataset ##

cow <- read.dta13("data/cow/NMC_5_0.dta")

summary(cow$upop)
summary(cow$tpop)

cow <- cow %>%
  mutate(upop = ifelse(upop <0, NA, upop), 
         urban_percent = upop/tpop) %>%
  rename(cown = ccode) %>%
  select(cown, year, urban_percent)

summary(cow$urban_percent)

#### Additional Data for Hong Kong #### 
world_bank_hong_kong <- read.csv2("data/world_bank_data/hong_kong_urban2.csv", na.strings = "..")

world_bank_hong_kong$Urban.population....of.total.population. <- as.numeric(as.character(world_bank_hong_kong$Urban.population....of.total.population.))

world_bank_hong_kong <- world_bank_hong_kong %>%
  rename(year = Time,
         country_name = Country.Name, 
         urban_percent = Urban.population....of.total.population.) %>%
  select(year, urban_percent) %>%
  mutate(urban_percent = urban_percent/100)

world_bank_hong_kong$cown <- 715

cow <- rbind(cow, world_bank_hong_kong)
cow$cown[cow$cown == 816] <- 817


## World Bank Data on Manufacturing Sector ##

# Manufacturing, value added (% of GDP)

worldbank_data <- read.csv("data/world_bank_data/730db413-06a3-4f58-810a-f9341fd46150_Data.csv", na.strings = "..")

worldbank_data <- worldbank_data %>%
  rename(year = �..Time, 
         country_name = Country.Name, 
         manufacturing_percent = Manufacturing..value.added....of.GDP...NV.IND.MANF.ZS.) %>%
  select(year, country_name, manufacturing_percent)

worldbank_data <- worldbank_data[1:15840,]

worldbank_data$year <- as.numeric(as.character(worldbank_data$year))

worldbank_data$cown <- countrycode(worldbank_data$country_name, "country.name", "cown", warn = TRUE)

worldbank_data$cown[worldbank_data$country_name == "Hong Kong SAR, China"] <- 715

worldbank_data <- worldbank_data %>%
  drop_na(cown) %>%
  select(-country_name)

#### Merge Data ####

vdem <- vdem %>%
  left_join(swiid_summary, c("cown", "year"))

vdem <- vdem %>%
  left_join(cow, c("cown", "year"))

vdem <- vdem %>%
  left_join(worldbank_data, c("cown", "year"))

summary(vdem$cown)

## Interpolate Missing Values of Manufacturing Percent GDP ##

vdem_man <- vdem %>%
  group_by(country_id) %>%
  filter(any(!is.na(manufacturing_percent))) %>%
  select(country_id, year, manufacturing_percent )

vdem_man <- vdem_man %>%
  group_by(country_id) %>%
  mutate(manufacturing_percent_ipol = na_interpolation(manufacturing_percent, option = "spline")) %>%
  select(-manufacturing_percent)

vdem <- vdem %>%
  left_join(vdem_man, c("country_id", "year"))


#### Building Variables political inclusiveness and economic inclusiveness ####

vdem <- vdem %>%
  mutate(pol_incl = (v2pepwrses+v2pepwrsoc)/2)
summary(vdem$pol_incl) #smaller values indicate more exclusionary regimes, higher values indicate more inclusionary regimes

vdem <- vdem %>%
  drop_na(v2x_regime)

#### Generate variables ####

vdem <- vdem %>%
  mutate(e_wb_pop_ln = log10(e_wb_pop),
         e_total_oil_income_pc = ifelse(e_total_oil_income_pc==0, 1, e_total_oil_income_pc),
         e_total_oil_income_pc_ln = log10(e_total_oil_income_pc))

#### Generate Competitive and Hegemonic Multiparty regimes ####

## prepare Election specific variables ##

summary(vdem$v2elmulpar_ord)
summary(vdem$v2elfrfair_ord)

vdem <- vdem %>%
  group_by(country_id) %>%
  fill(v2elmulpar_ord) %>%
  fill(v2elfrfair_ord)

vdem <- vdem %>%
  mutate(v2elfrfair_ord = if_else(is.na(v2elfrfair_ord), 0, v2elfrfair_ord), # replace NA with = 0, assumption: NA no elections in place
         v2elmulpar_ord = if_else(is.na(v2elmulpar_ord), 0, v2elmulpar_ord))


summary(vdem$v2elmulpar_ord)
summary(vdem$v2elfrfair_ord)
summary(vdem$v2elsuffrage)


vdem <- vdem %>%
  mutate(auto_regime_type = case_when(v2elmulpar_ord >=2 & v2elfrfair_ord >= 2 & v2elsuffrage > 25 ~ 2, # competitive mulitparty regime
                                      v2elmulpar_ord >=2 & v2elfrfair_ord < 2 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord >=2 & v2elsuffrage <= 25 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord <2  ~ 0 )) # closed autocracy

table(vdem$auto_regime_type) 
summary(vdem$auto_regime_type)


#### Generate Sample ####

vdem <- vdem %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_rel <- vdem %>%
  drop_na(rel_red)

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE) # Distinct observations

vdem_rel <- vdem_rel %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_rel <- vdem_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_rel2 <- vdem_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_rel <- vdem_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, rel_red, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_rel2 <- vdem_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, rel_red, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

#### Computing group_mean and de-meaned variables ####

vdem_rel <- sjmisc::de_mean(vdem_rel,pol_incl, v2xps_party, v2x_regime,
                            v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                            grp = country_id)

vdem_rel2 <- sjmisc::de_mean(vdem_rel2, pol_incl, v2xps_party, v2x_regime,
                             e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

#### Descriptive statistics, Appendix A1 ##

cols1 <- c("rel_red", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols2 <- c("rel_red", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 


## reference Level for Factor ##
vdem_rel$auto_regime_type  <- as.factor(vdem_rel$auto_regime_type)
vdem_rel2$auto_regime_type  <- as.factor(vdem_rel2$auto_regime_type)

vdem_rel$auto_regime_type = relevel(vdem_rel$auto_regime_type, ref=1)
vdem_rel2auto_regime_type = relevel(vdem_rel2$auto_regime_type, ref=1)


# Model 1 #
stargazer(
  title="Summary Statistics for Redistribution Dataset", 
  as.data.frame(vdem_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Rel. Redistribution", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 2-3 #
stargazer(
  title="Summary Statistics for Redistribution Dataset", 
  as.data.frame(vdem_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Rel. Redistribution", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

#### With binary regime indicator ####

m1 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE

m2 <- lmer(rel_red ~  year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
             v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
             auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
             v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
             v2pscohesv_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_rel, 
           REML = TRUE, 
           control=lmerControl(optCtrl=list(maxfun=1e9)))
isSingular(m2, tol = 1e-05) # FALSE

m3 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m3, tol = 1e-05) # FALSE

m4 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
             v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
             auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
             v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
             v2pscohesv_dm*auto_regime_type + 
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m4, tol = 1e-05) # FALSE

# Interaction Models, compare Giesselmann and Schmidt Catran 2018 (double-demeanining of interactions)

m5 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1 | country_id),
           data = vdem_rel2)
isSingular(m5, tol = 1e-05) # FALSE

m6 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
             v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
             auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
             v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
             v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m6, tol = 1e-05) # FALSE

tab_model(m1, m2, m3, m4, m5,  
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"))

#### Regression Output ####

texreg(list(m1, m3, m5), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "tab:tab-01", 
       longtable = TRUE)

texreg(list(m2, m4, m6), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party organizations", 
                             "Party organizations (b)", 
                             "Party Branches", 
                             "Party Branches (b)", 
                             "Party Linkages", 
                             "Party Linkages (b)", 
                             "Distinct Party Platforms", 
                             "Distinct Party Platforms (b)", 
                             "Legislative Party Cohesion", 
                             "Legislative Party Cohesion (b)", 
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * Party organizations", 
                             "CM Autocracy * Party organizations", 
                             "HM Autocracy * Party Branches", 
                             "CM Autocracy * Party Branches", 
                             "HM Autocracy * Party Linkages", 
                             "CM Autocracy * Party Linkages", 
                             "HM Autocracy * Distinct Party Platforms", 
                             "CM Autocracy * Distinct Party Platforms", 
                             "HM Autocracy * Legislative Party Cohesion",
                             "CM Autocracy * Legislative Party Cohesion",
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "Inclusiveness * Party organizations",
                             "Inclusiveness * Party Branches",
                             "Inclusiveness * Party Linkages",
                             "Inclusiveness * Distinct Party Platforms",
                             "Inclusiveness * Legislative Party Cohesion",
                             "Inclusiveness * Party organizations * HM Autocracy",
                             "Inclusiveness * Party organizations * CM Autocracy",
                             "Inclusiveness * Party Branches * HM Autocracy",
                             "Inclusiveness * Party Branches * CM Autocracy",
                             "Inclusiveness * Party Linkages * HM Autocracy",
                             "Inclusiveness * Party Linkages * CM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * HM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * CM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * HM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * CM Autocracy"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Redistribution / Universalism",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:B2", 
       longtable = TRUE)

#### Marginal Plots for m3, m4 and m5 and m6 ####


#### Figure 1  based on Model 3 ####

#two way intercation effect plots 

margin_data <- ggpredict(m3, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)

ggsave("Output/Margins/S1_Figure_1.pdf", height = 20, width = 24, units= c("cm"))

#### Figure 2 Interaction effects ####

sd(vdem_rel2$pol_incl_dm)

margin_data <- ggpredict(m5, terms = c("v2xps_party_dm", "pol_incl_dm[-0.33, 0.33]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.33", "0.33"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)


ggsave("Output/Margins/S1_Figure_2.pdf", height = 12, width = 24, units= c("cm"))


##########################################################################
##########################################################################

#### APPENDIX c ####
#### 2. REWB with democratic and autocratic countries  ####

# clear workspace
rm(list=ls())

#### Load Data ####
vdem <- readRDS("data/vdem_merged.rds") 
vdem_spaw <- readRDS("data/vdem_merged_spaw.rds") 

vdem <- vdem %>%
  drop_na(v2x_regime)

vdem_spaw <- vdem_spaw %>%
  drop_na(v2x_regime)

#### Generate Regime Types ####

vdem <- vdem %>%
  mutate(auto_regime_type = case_when(v2x_regime >=2 ~ 3, # liberal and electoral democracies
                                      v2elmulpar_ord >=2 & v2elfrfair_ord >= 2 & v2elsuffrage > 25 ~ 2, # competitive mulitparty regime
                                      v2elmulpar_ord >=2 & v2elfrfair_ord < 2 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord >=2 & v2elsuffrage <= 25 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord <2  ~ 0 )) # closed autocracy

table(vdem$auto_regime_type) 
summary(vdem$auto_regime_type)


vdem_spaw <- vdem_spaw %>%
  mutate(auto_regime_type = case_when(v2x_regime >=2 ~ 3, # liberal and electoral democracies
                                      v2elmulpar_ord >=2 & v2elfrfair_ord >= 2 & v2elsuffrage > 25 ~ 2, # competitive mulitparty regime
                                      v2elmulpar_ord >=2 & v2elfrfair_ord < 2 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord >=2 & v2elsuffrage <= 25 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord <2  ~ 0 )) # closed autocracy

table(vdem_spaw$auto_regime_type) 
summary(vdem_spaw$auto_regime_type)

#### Generate Sample ####
vdem_rel <- vdem %>%
  drop_na(rel_red)

vdem_spaw <- vdem_spaw %>%
  drop_na(universalism_all)

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE) # Distinct observations
vdem_spaw <- distinct(vdem_spaw, country_id, year, .keep_all= TRUE) # Distinct observations

vdem_rel <- vdem_rel %>% # delete those countries with less than 5 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_spaw <- vdem_spaw %>% # delete those countries with less than 5 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)
vdem_rel <- vdem_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_spaw <- vdem_spaw %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_rel2 <- vdem_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_spaw2 <- vdem_spaw %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_rel <- vdem_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, rel_red, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_rel2 <- vdem_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, rel_red, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

vdem_spaw <- vdem_spaw %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, universalism_all, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_spaw2 <- vdem_spaw2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, universalism_all, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

#### Computing group_mean and de-meaned variables ####

vdem_rel <- sjmisc::de_mean(vdem_rel,pol_incl, v2xps_party, v2x_regime,
                            v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                            grp = country_id)

vdem_rel2 <- sjmisc::de_mean(vdem_rel2, pol_incl, v2xps_party, v2x_regime,
                             e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

vdem_spaw <- sjmisc::de_mean(vdem_spaw,pol_incl, v2xps_party, v2x_regime,
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

vdem_spaw2 <- sjmisc::de_mean(vdem_spaw2, pol_incl, v2xps_party, v2x_regime,
                              e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                              v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                              grp = country_id)

#### Descriptive statistics, Appendix C1 ##

cols1 <- c("rel_red", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols2 <- c("rel_red", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols3 <- c("universalism_all", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols4 <- c("universalism_all", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm")

# Model 1 #
stargazer(
  title="Summary Statistics for Redistribution Dataset", 
  as.data.frame(vdem_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Rel. Redistribution", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 2 #
stargazer(
  title="Summary Statistics for Universalism Dataset", 
  as.data.frame(vdem_spaw[, cols3]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Universalism", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 3 and 5 #
stargazer(
  title="Summary Statistics for Redistribution Dataset", 
  as.data.frame(vdem_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Rel. Redistribution", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 4 and 6 #
stargazer(
  title="Summary Statistics for Universalism Dataset", 
  as.data.frame(vdem_spaw2[, cols4]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Universalism", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)


## reference Level for Factor ##
vdem_rel$auto_regime_type  <- as.factor(vdem_rel$auto_regime_type)
vdem_rel2$auto_regime_type  <- as.factor(vdem_rel2$auto_regime_type)
vdem_rel$auto_regime_type = relevel(vdem_rel$auto_regime_type, ref=1)
vdem_rel2$auto_regime_type = relevel(vdem_rel2$auto_regime_type, ref=1)

vdem_spaw$auto_regime_type  <- as.factor(vdem_spaw$auto_regime_type)
vdem_spaw2$auto_regime_type  <- as.factor(vdem_spaw2$auto_regime_type)
vdem_spaw$auto_regime_type = relevel(vdem_spaw$auto_regime_type, ref=1)
vdem_spaw2$auto_regime_type = relevel(vdem_spaw2$auto_regime_type, ref=1)

#### Models ####
#### Main Models ####

m1 <- lmer(rel_red ~  year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE
performance::check_model(m1)


m2 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_spaw, 
           REML = TRUE) # REML estimation for parameters
isSingular(m2, tol = 1e-05) # FALSE
performance::check_model(m2)

m3 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m3, tol = 1e-05) # FALSE
performance::check_model(m3)


m4 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_spaw2)
isSingular(m4, tol = 1e-05) # FALSE
performance::check_model(m4)


m5 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m5, tol = 1e-05) # FALSE
performance::check_model(m5)


m6 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_spaw2)
isSingular(m6, tol = 1e-05) # FALSE
performance::check_model(m6)


tab_model(m1, m2, m3, m4, m5, m6, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"))

#### Regression Output ####

texreg(list(m1, m2, m3, m4, m5, m6), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2, custom.coef.names = c("(Intercept)",
                                                   "Year",
                                                   "Pol Inclusiveness",
                                                   "Pol Inclusiveness (b)",
                                                   "Party Institutionalization",
                                                   "Party Institutionalization (b)",
                                                   "Hegemonic Multiparty Autocracy",
                                                   "Competitive Multiparty Autocracy",
                                                   "Democracy", 
                                                   "HM Autocracy * Inclusiveness",
                                                   "CM Autocracy * Inclusiveness",
                                                   "Democracy * Inclusiveness",
                                                   "HM Autocracy * PI",
                                                   "CM Autocracy  * PI", 
                                                   "Democracy  * PI", 
                                                   "GDP pc log",
                                                   "GDP pc log (b)",
                                                   "Population log",
                                                   "Population log (b)",
                                                   "Communist Ideology", 
                                                   "Communist Ideology (b)", 
                                                   "Urban Pop %",
                                                   "Urban Pop % (b)",
                                                   "Manufacturing Sector %",
                                                   "Manufacturing Sector % (b)",
                                                   "PI * Inclusiveness",
                                                   "HM Autocracy * PI * Inclusiveness",
                                                   "CM Autocracy * PI * Inclusiveness", 
                                                   "Democracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:C5", 
       longtable = TRUE)

#### Marginal Plots for m3, m4 and m5 and m6 ####
#### Figure 1  based on Model 3 and 4 ####
#two way intercation effect plots 

margin_data <- ggpredict(m3, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2", "3"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy",
                                       "Competitive Multiparty Autocracy", "Democracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2", "3"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy",
                                       "Competitive Multiparty Autocracy", "Democracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m4, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2", "3"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy",
                                       "Competitive Multiparty Autocracy", "Democracy"))

plot3 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "Predicted universalism old-age pensions") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m4, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2", "3"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy",
                                       "Competitive Multiparty Autocracy", "Democracy"))

plot4 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted universalism old-age pensions") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


ggarrange(plot1, plot2, plot3, plot4,
          labels = c("a", "b", "c", "d"), 
          ncol = 2, nrow = 2, 
          common.legend = TRUE)

ggsave("Output/Margins/S2_Figure_1.pdf", height =20, width = 24, units= c("cm"))

#### Figure 2 Three-Way Interaction effects ####

sd(vdem_rel2$pol_incl_dm)

margin_data <- ggpredict(m5, terms = c("v2xps_party_dm", "pol_incl_dm[-0.39, 0.39]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2", "3"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy",
                                       "Competitive Multiparty Autocracy", "Democracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.39", "0.39"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)



sd(vdem_spaw2$pol_incl_dm)

margin_data <- ggpredict(m6, terms = c("v2xps_party_dm", "pol_incl_dm[-0.43, 0.43]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2", "3"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy",
                                       "Competitive Multiparty Autocracy", "Democracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.43", "0.43"),
                            labels = c("-1 SD", "+ 1 SD"))

plot2 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted universalism old-age pensions") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)


ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)

ggsave("Output/Margins/S2_Figure_2.pdf", height = 25, width = 20, units= c("cm"))

##########################################################################
##########################################################################

#### APPENDIX D ####
#### 3. Albertus and Menaldo (2014) data  ####

rm(list=ls())


## V-DEM
vdem <- read_csv("data/vdem_9/V-Dem-CY-Full+Others-v9.csv") 

vdem$cown <- countrycode(vdem$country_name, "country.name", "cown", warn = TRUE)
vdem$cown[vdem$country_name == "Hong Kong"] <- 715
vdem$cown[vdem$country_name == "Republic of Vietnam"] <- 817
vdem$cown <- as.integer(vdem$cown)


sp_data <- read_csv("data/sp_data/AMGaming_DATA.csv") 
sp_data$cown <- countrycode(sp_data$cnamehabmen, "country.name", "cown", warn = TRUE)
sp_data$cown[sp_data$cnamehabmen == "Hong Kong"] <- "715"
sp_data$cown <- as.integer(sp_data$cown)


## Correlates of War National Material Capabilities Dataset ##

cow <- read.dta13("data/cow/NMC_5_0.dta")

summary(cow$upop)
summary(cow$tpop)

cow <- cow %>%
  mutate(upop = ifelse(upop <0, NA, upop), 
         urban_percent = upop/tpop) %>%
  rename(cown = ccode) %>%
  select(cown, year, urban_percent)

summary(cow$urban_percent)

#### Additional Data for Hong Kong #### 
world_bank_hong_kong <- read.csv2("data/world_bank_data/hong_kong_urban2.csv", na.strings = "..")

world_bank_hong_kong$Urban.population....of.total.population. <- as.numeric(as.character(world_bank_hong_kong$Urban.population....of.total.population.))

world_bank_hong_kong <- world_bank_hong_kong %>%
  rename(year = Time,
         country_name = Country.Name, 
         urban_percent = Urban.population....of.total.population.) %>%
  select(year, urban_percent) %>%
  mutate(urban_percent = urban_percent/100)

world_bank_hong_kong$cown <- 715

cow <- rbind(cow, world_bank_hong_kong)
cow$cown[cow$cown == 816] <- 817

## World Bank Data on Manufacturing Sector ##

# Manufacturing, value added (% of GDP)

worldbank_data <- read.csv("data/world_bank_data/730db413-06a3-4f58-810a-f9341fd46150_Data.csv", na.strings = "..")

worldbank_data <- worldbank_data %>%
  rename(year = �..Time, 
         country_name = Country.Name, 
         manufacturing_percent = Manufacturing..value.added....of.GDP...NV.IND.MANF.ZS.) %>%
  select(year, country_name, manufacturing_percent)

worldbank_data <- worldbank_data[1:15840,]

worldbank_data$year <- as.numeric(as.character(worldbank_data$year))

worldbank_data$cown <- countrycode(worldbank_data$country_name, "country.name", "cown", warn = TRUE)

worldbank_data$cown[worldbank_data$country_name == "Hong Kong SAR, China"] <- 715

worldbank_data <- worldbank_data %>%
  drop_na(cown) %>%
  select(-country_name)

#### Merge Datasets ####

vdem <- vdem %>%
  left_join(sp_data, c("cown", "year"))

vdem <- vdem %>%
  left_join(cow, c("cown", "year"))

vdem <- vdem %>%
  left_join(worldbank_data, c("cown", "year"))

## Interpolate Missing Values of Manufacturing Percent GDP ##

vdem_man <- vdem %>%
  group_by(country_id) %>%
  filter(any(!is.na(manufacturing_percent))) %>%
  select(country_id, year, manufacturing_percent )

vdem_man <- vdem_man %>%
  group_by(country_id) %>%
  mutate(manufacturing_percent_ipol = na_interpolation(manufacturing_percent, option = "spline")) %>%
  select(-manufacturing_percent)

vdem <- vdem %>%
  left_join(vdem_man, c("country_id", "year"))

#### Building Variables political inclusiveness and economic inclusiveness ####

vdem <- vdem %>%
  mutate(pol_incl = (v2pepwrses+v2pepwrsoc)/2,
         eco_incl = v2dlencmps )
summary(vdem$pol_incl) #smaller values indicate more exclusionary regimes, higher values indicate more inclusionary regimes
summary(vdem$eco_incl) #smaller values indicate more exclusionary regimes, higher values indicate more inclusionary regimes

vdem <- vdem %>%
  mutate(incl = (pol_incl + eco_incl)/2) #smaller values indicate more exclusionary regimes, higher values indicate more inclusionary regimes
summary(vdem$incl) 

#### Generate variables ####

vdem <- vdem %>%
  mutate(e_wb_pop_ln = log10(e_wb_pop),
         e_total_oil_income_pc = ifelse(e_total_oil_income_pc==0, 1, e_total_oil_income_pc),
         e_total_oil_income_pc_ln = log10(e_total_oil_income_pc))

## Percent Urban Population ##

vdem <- vdem %>%
  group_by(country_id) %>%
  fill(urban_percent, .direction = "updown") # Fill missing oberservations with next or previos value. Affect 157 NA in 2013-2018

vdem$manufacturing_percent_ipol <- vdem$manufacturing_percent_ipol/100

#### Generate Competitive and Hegemonic Multiparty regimes ####

## prepare Election specific variables ##

summary(vdem$v2elmulpar_ord)
summary(vdem$v2elfrfair_ord)

vdem <- vdem %>%
  group_by(country_id) %>%
  fill(v2elmulpar_ord) %>%
  fill(v2elfrfair_ord)

vdem <- vdem %>%
  mutate(v2elfrfair_ord = if_else(is.na(v2elfrfair_ord), 0, v2elfrfair_ord), # replace NA with = 0, assumption: NA no elections in place
         v2elmulpar_ord = if_else(is.na(v2elmulpar_ord), 0, v2elmulpar_ord))


summary(vdem$v2elmulpar_ord)
summary(vdem$v2elfrfair_ord)
summary(vdem$v2elsuffrage)


vdem <- vdem %>%
  mutate(auto_regime_type = case_when(v2elmulpar_ord >=2 & v2elfrfair_ord >= 2 & v2elsuffrage > 25 ~ 2, # competitive mulitparty regime
                                      v2elmulpar_ord >=2 & v2elfrfair_ord < 2 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord >=2 & v2elsuffrage <= 25 ~ 1, # hegemonic mulitparty regime
                                      v2elmulpar_ord <2  ~ 0 )) # closed autocracy

table(vdem$auto_regime_type) 
summary(vdem$auto_regime_type)

#### Export V-Dem-SWIID DATASET ####

# export basic dataset 
saveRDS(vdem, file = "data/vdem_social.rds")
vdem_social <- readRDS("data/vdem_social.rds")

summary(vdem_social$inc_gdp_ratio) # Taxes on Income and profits and capital gains in % GDP
summary(vdem_social$social_spending_gdp_perc) # Social spending in % GDP
summary(vdem_social$social_protection_pct) # Welfare and Insurance spending in % GDP


########################################################
#### Income Taxation in % GDP, APPENDIX D1
########################################################

#### Cleaning data and observations ####

vdem_social <- vdem_social %>%
  drop_na(v2x_regime)

vdem_social <- vdem_social %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_social_rel <- vdem_social %>%
  drop_na(inc_gdp_ratio) ###  social spending, change for income taxation inc_gdp_ratio

vdem_social_rel <- distinct(vdem_social_rel, country_id, year, .keep_all= TRUE)

vdem_social_rel <- vdem_social_rel %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_social_rel <- vdem_social_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_social_rel2 <- vdem_social_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1 )

vdem_social_rel <- vdem_social_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, inc_gdp_ratio, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_social_rel2 <- vdem_social_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, inc_gdp_ratio, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)


#### Computing group_mean and de-meaned variables ####

vdem_social_rel <- sjmisc::de_mean(vdem_social_rel, pol_incl, v2xps_party, v2x_regime,
                                   v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                                   grp = country_id)

vdem_social_rel2 <- sjmisc::de_mean(vdem_social_rel2, pol_incl, v2xps_party, v2x_regime,
                                    e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                                    v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                                    grp = country_id)

#### Descriptive statistics, Appendix A1 ##

cols1 <- c("inc_gdp_ratio", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 


cols2 <- c("inc_gdp_ratio", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

# Model 1 and 4 #
stargazer(
  title="Summary Statistics for Taxes on Income, Profits, and Capital Gains in % of GDP Dataset", 
  as.data.frame(vdem_social_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Income GDP ratio", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)


# Model 2-3, 5-6 #

stargazer(
  title="Summary Statistics for Taxes on Income, Profits, and Capital Gains in % of GDP Dataset", 
  as.data.frame(vdem_social_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Income GDP ratio", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion") )


## reference Level for Factor ##
vdem_social_rel$auto_regime_type  <- as.factor(vdem_social_rel$auto_regime_type)
vdem_social_rel2$auto_regime_type  <- as.factor(vdem_social_rel2$auto_regime_type)

vdem_social_rel$auto_regime_type = relevel(vdem_social_rel$auto_regime_type, ref=1)
vdem_social_rel2$auto_regime_type = relevel(vdem_social_rel2$auto_regime_type, ref=1)


#### Models ####

m1 <- lmer(inc_gdp_ratio ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_social_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE

m1.1 <- lmer(inc_gdp_ratio ~  year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               (1  | country_id),
           data = vdem_social_rel, 
           REML = TRUE, 
           control=lmerControl(optCtrl=list(maxfun=1e9)))
isSingular(m1.1, tol = 1e-05) # FALSE

m2 <- lmer(inc_gdp_ratio ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_social_rel2)
isSingular(m2, tol = 1e-05) # FALSE

m2.1 <- lmer(inc_gdp_ratio ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
           data = vdem_social_rel2)
isSingular(m2.1, tol = 1e-05) # FALSE

# Interaction Models, compare Giesselmann and Schmidt Catran 2018 (double-demeanining of interactions)

m3 <- lmer(inc_gdp_ratio ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_social_rel2)
isSingular(m3, tol = 1e-05) # FALSE

m3.1 <- lmer(inc_gdp_ratio ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
               v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
               v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
           data = vdem_social_rel2)
isSingular(m3.1, tol = 1e-05) # FALSE

tab_model(m1, m2, m3, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3"))

#### Regression Output ####

texreg(list(m1, m2, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Income Taxation in % GDP",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:D3", 
       longtable = TRUE)


texreg(list(m1.1, m2.1, m3.1), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party organizations", 
                             "Party organizations (b)", 
                             "Party Branches", 
                             "Party Branches (b)", 
                             "Party Linkages", 
                             "Party Linkages (b)", 
                             "Distinct Party Platforms", 
                             "Distinct Party Platforms (b)", 
                             "Legislative Party Cohesion", 
                             "Legislative Party Cohesion (b)", 
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * Party organizations", 
                             "CM Autocracy * Party organizations", 
                             "HM Autocracy * Party Branches", 
                             "CM Autocracy * Party Branches", 
                             "HM Autocracy * Party Linkages", 
                             "CM Autocracy * Party Linkages", 
                             "HM Autocracy * Distinct Party Platforms", 
                             "CM Autocracy * Distinct Party Platforms", 
                             "HM Autocracy * Legislative Party Cohesion",
                             "CM Autocracy * Legislative Party Cohesion",
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "Inclusiveness * Party organizations",
                             "Inclusiveness * Party Branches",
                             "Inclusiveness * Party Linkages",
                             "Inclusiveness * Distinct Party Platforms",
                             "Inclusiveness * Legislative Party Cohesion",
                             "Inclusiveness * Party organizations * HM Autocracy",
                             "Inclusiveness * Party organizations * CM Autocracy",
                             "Inclusiveness * Party Branches * HM Autocracy",
                             "Inclusiveness * Party Branches * CM Autocracy",
                             "Inclusiveness * Party Linkages * HM Autocracy",
                             "Inclusiveness * Party Linkages * CM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * HM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * CM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * HM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * CM Autocracy"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Income Taxation in % GDP",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:D4", 
       longtable = TRUE)


#### Marginal Plots ####

#### Figure 1, based on Models 3 ####

margin_data <- ggpredict(m2, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "redicted Income Taxation") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_social_rel2$inc_gdp_ratio), linetype="dashed")


margin_data <- ggpredict(m2, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted Income Taxation") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_social_rel2$inc_gdp_ratio), linetype="dashed")


ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)

ggsave("output/Margins/S3_Figure_1.pdf", height = 20, width = 24, units= c("cm"))

#### Figure 2 Three-Way Interaction effects ####

sd(vdem_social_rel2$pol_incl_dm)

margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "pol_incl_dm[-0.37, 0.37]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.37", "0.37"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted Income Taxation") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_social_rel2$inc_gdp_ratio), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)


ggarrange(plot1, 
          labels = c("a"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)
ggsave("output/Margins/S3_Figure_2.pdf", height = 14, width = 24, units= c("cm"))

########################################################
#### 4. Social Spending Data APPENDIX D2
########################################################

#### Cleaning data and observations ####

vdem_social <- readRDS("data/vdem_social.rds")

vdem_social <- vdem_social %>%
  drop_na(v2x_regime)

vdem_social <- vdem_social %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_social_rel <- vdem_social %>%
  drop_na(social_spending_gdp_perc) ## social spending

vdem_social_rel <- distinct(vdem_social_rel, country_id, year, .keep_all= TRUE)

vdem_social_rel <- vdem_social_rel %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_social_rel <- vdem_social_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_social_rel2 <- vdem_social_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol )

vdem_social_rel <- vdem_social_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, social_spending_gdp_perc, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_social_rel2 <- vdem_social_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, social_spending_gdp_perc, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)


#### Computing group_mean and de-meaned variables ####

vdem_social_rel <- sjmisc::de_mean(vdem_social_rel, pol_incl, v2xps_party, v2x_regime,
                                   v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                                   grp = country_id)

vdem_social_rel2 <- sjmisc::de_mean(vdem_social_rel2, pol_incl, v2xps_party, v2x_regime,
                                    e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                                    v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                                    grp = country_id)


#### Descriptive statistics, Appendix D2 ##

cols1 <- c("social_spending_gdp_perc",  "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols2 <- c("social_spending_gdp_perc", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm")

# Model 1 and 1.1 #
stargazer(
  title="Summary Statistics for Social Spending Dataset", 
  as.data.frame(vdem_social_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Social Spending % GDP", 
                      "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 2-3 and 2.1 and 3.1 #

stargazer(
  title="Summary Statistics for Social Spending Dataset", 
  as.data.frame(vdem_social_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Social Spending % GDP", 
                      "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

## reference Level for Factor ##
vdem_social_rel$auto_regime_type  <- as.factor(vdem_social_rel$auto_regime_type)
vdem_social_rel2$auto_regime_type  <- as.factor(vdem_social_rel2$auto_regime_type)

vdem_social_rel$auto_regime_type = relevel(vdem_social_rel$auto_regime_type, ref=1)
vdem_social_rel2$auto_regime_type = relevel(vdem_social_rel2$auto_regime_type, ref=1)


#### Models ####

m1 <- lmer(social_spending_gdp_perc ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_social_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE

m1.1 <- lmer(social_spending_gdp_perc ~  year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               (1  | country_id),
             data = vdem_social_rel, 
             REML = TRUE, 
             control=lmerControl(optCtrl=list(maxfun=1e9)))
isSingular(m1.1, tol = 1e-05) # FALSE

m2 <- lmer(social_spending_gdp_perc ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_social_rel2)
isSingular(m2, tol = 1e-05) # FALSE

m2.1 <- lmer(social_spending_gdp_perc ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
             data = vdem_social_rel2)
isSingular(m2.1, tol = 1e-05) # FALSE

# Interaction Models, compare Giesselmann and Schmidt Catran 2018 (double-demeanining of interactions)

m3 <- lmer(social_spending_gdp_perc ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_social_rel2)
isSingular(m3, tol = 1e-05) # FALSE

m3.1 <- lmer(social_spending_gdp_perc ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
               v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
               v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
             data = vdem_social_rel2)
isSingular(m3.1, tol = 1e-05) # FALSE

tab_model(m1, m2, m3, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3"))

#### Regression Output ####

texreg(list(m1, m2, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Social Spending per Capita in % GDP",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:D7", 
       longtable = TRUE)


texreg(list(m1.1, m2.1, m3.1), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party organizations", 
                             "Party organizations (b)", 
                             "Party Branches", 
                             "Party Branches (b)", 
                             "Party Linkages", 
                             "Party Linkages (b)", 
                             "Distinct Party Platforms", 
                             "Distinct Party Platforms (b)", 
                             "Legislative Party Cohesion", 
                             "Legislative Party Cohesion (b)", 
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * Party organizations", 
                             "CM Autocracy * Party organizations", 
                             "HM Autocracy * Party Branches", 
                             "CM Autocracy * Party Branches", 
                             "HM Autocracy * Party Linkages", 
                             "CM Autocracy * Party Linkages", 
                             "HM Autocracy * Distinct Party Platforms", 
                             "CM Autocracy * Distinct Party Platforms", 
                             "HM Autocracy * Legislative Party Cohesion",
                             "CM Autocracy * Legislative Party Cohesion",
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "Inclusiveness * Party organizations",
                             "Inclusiveness * Party Branches",
                             "Inclusiveness * Party Linkages",
                             "Inclusiveness * Distinct Party Platforms",
                             "Inclusiveness * Legislative Party Cohesion",
                             "Inclusiveness * Party organizations * HM Autocracy",
                             "Inclusiveness * Party organizations * CM Autocracy",
                             "Inclusiveness * Party Branches * HM Autocracy",
                             "Inclusiveness * Party Branches * CM Autocracy",
                             "Inclusiveness * Party Linkages * HM Autocracy",
                             "Inclusiveness * Party Linkages * CM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * HM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * CM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * HM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * CM Autocracy"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Social Spending per Capita in % GDP",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:D8", 
       longtable = TRUE)

#### Marginal Plots for m3, m4 and m5 and m6 ####


#### Figure D3 ####

margin_data <- ggpredict(m2, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "redicted Social Spending") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_social_rel2$social_spending_gdp_perc), linetype="dashed")


margin_data <- ggpredict(m2, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted Social Spending") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_social_rel2$social_spending_gdp_perc), linetype="dashed")


ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)
ggsave("output/Margins/S4_Figure_1.pdf", height = 20, width = 24, units= c("cm"))

#### Figure D4 Three-Way Interaction effects ####

sd(vdem_social_rel2$pol_incl_dm)

margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "pol_incl_dm[-0.33, 0.33]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.33", "0.33"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted Social Spending") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_social_rel2$social_spending_gdp_perc), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)

ggarrange(plot1, 
          ncol = 1, 
          common.legend = TRUE)
ggsave("output/Margins/S4_Figure_2.pdf", height = 16, width = 24, units= c("cm"))

########################################################
#### 4. Social Protection Data APPENDIX D3
########################################################

#### Cleaning data and observations ####

vdem_social <- readRDS("data/vdem_social.rds")

vdem_social <- vdem_social %>%
  drop_na(v2x_regime)

vdem_social <- vdem_social %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_social_rel <- vdem_social %>%
  drop_na(social_protection_pct) ## social spending

vdem_social_rel <- distinct(vdem_social_rel, country_id, year, .keep_all= TRUE)

vdem_social_rel <- vdem_social_rel %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_social_rel <- vdem_social_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_social_rel2 <- vdem_social_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol )

vdem_social_rel <- vdem_social_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, social_protection_pct, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_social_rel2 <- vdem_social_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, social_protection_pct, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)


#### Computing group_mean and de-meaned variables ####

vdem_social_rel <- sjmisc::de_mean(vdem_social_rel, pol_incl, v2xps_party, v2x_regime,
                                   v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                                   grp = country_id)

vdem_social_rel2 <- sjmisc::de_mean(vdem_social_rel2, pol_incl, v2xps_party, v2x_regime,
                                    e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                                    v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                                    grp = country_id)


#### Descriptive statistics, Appendix D2 ##

cols1 <- c("social_protection_pct",  "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols2 <- c("social_protection_pct", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm")

# Model 1 and 1.1 #
stargazer(
  title="Summary Statistics for Social Protection Dataset", 
  as.data.frame(vdem_social_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Social Spending % GDP", 
                      "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 2-3 and 2.1 and 3.1 #

stargazer(
  title="Summary Statistics for Social Protection Dataset", 
  as.data.frame(vdem_social_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Social Spending % GDP", 
                      "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

## reference Level for Factor ##
vdem_social_rel$auto_regime_type  <- as.factor(vdem_social_rel$auto_regime_type)
vdem_social_rel2$auto_regime_type  <- as.factor(vdem_social_rel2$auto_regime_type)

vdem_social_rel$auto_regime_type = relevel(vdem_social_rel$auto_regime_type, ref=1)
vdem_social_rel2$auto_regime_type = relevel(vdem_social_rel2$auto_regime_type, ref=1)


#### Models ####

m1 <- lmer(social_protection_pct ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_social_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE

m1.1 <- lmer(social_protection_pct ~  year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               (1  | country_id),
             data = vdem_social_rel, 
             REML = TRUE, 
             control=lmerControl(optCtrl=list(maxfun=1e9)))
isSingular(m1.1, tol = 1e-05) # FALSE

m2 <- lmer(social_protection_pct ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_social_rel2)
isSingular(m2, tol = 1e-05) # FALSE

m2.1 <- lmer(social_protection_pct ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
             data = vdem_social_rel2)
isSingular(m2.1, tol = 1e-05) # FALSE

# Interaction Models, compare Giesselmann and Schmidt Catran 2018 (double-demeanining of interactions)

m3 <- lmer(social_protection_pct ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_social_rel2)
isSingular(m3, tol = 1e-05) # FALSE

m3.1 <- lmer(social_protection_pct ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
               v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
               v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
             data = vdem_social_rel2)
isSingular(m3.1, tol = 1e-05) # FALSE

tab_model(m1, m2, m3, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3"))

#### Regression Output ####

texreg(list(m1, m2, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Welfare and Social Insurance Spending per Capita of % GDP",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:D11", 
       longtable = TRUE)


texreg(list(m1.1, m2.1, m3.1), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party organizations", 
                             "Party organizations (b)", 
                             "Party Branches", 
                             "Party Branches (b)", 
                             "Party Linkages", 
                             "Party Linkages (b)", 
                             "Distinct Party Platforms", 
                             "Distinct Party Platforms (b)", 
                             "Legislative Party Cohesion", 
                             "Legislative Party Cohesion (b)", 
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * Party organizations", 
                             "CM Autocracy * Party organizations", 
                             "HM Autocracy * Party Branches", 
                             "CM Autocracy * Party Branches", 
                             "HM Autocracy * Party Linkages", 
                             "CM Autocracy * Party Linkages", 
                             "HM Autocracy * Distinct Party Platforms", 
                             "CM Autocracy * Distinct Party Platforms", 
                             "HM Autocracy * Legislative Party Cohesion",
                             "CM Autocracy * Legislative Party Cohesion",
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "Inclusiveness * Party organizations",
                             "Inclusiveness * Party Branches",
                             "Inclusiveness * Party Linkages",
                             "Inclusiveness * Distinct Party Platforms",
                             "Inclusiveness * Legislative Party Cohesion",
                             "Inclusiveness * Party organizations * HM Autocracy",
                             "Inclusiveness * Party organizations * CM Autocracy",
                             "Inclusiveness * Party Branches * HM Autocracy",
                             "Inclusiveness * Party Branches * CM Autocracy",
                             "Inclusiveness * Party Linkages * HM Autocracy",
                             "Inclusiveness * Party Linkages * CM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * HM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * CM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * HM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * CM Autocracy"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Welfare and Social Insurance Spending per Capita of % GDP",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:D12", 
       longtable = TRUE)

#### Marginal Plots ###

#### Figure D5 ####

margin_data <- ggpredict(m2, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "redicted Social Protection") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_social_rel2$social_protection_pct), linetype="dashed")


margin_data <- ggpredict(m2, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted Social Protection") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_social_rel2$social_protection_pct), linetype="dashed")


ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)
ggsave("output/Margins/S5_Figure_1.pdf", height = 20, width = 24, units= c("cm"))

#### Figure 4 Interaction effects ####

sd(vdem_social_rel2$pol_incl_dm)

margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "pol_incl_dm[-0.33, 0.33]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.33", "0.33"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted Social Protection") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_social_rel2$social_protection_pct), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)

ggarrange(plot1, 
          ncol = 1, 
          common.legend = TRUE)
ggsave("output/Margins/S5_Figure_2.pdf", height = 14, width = 24, units= c("cm"))


##########################################################################
##########################################################################

#### APPENDIX E ####
#### 6. Public Goods Provision V-Dem Data ####

# clear workspace
rm(list=ls())

# set working directory

#### Load Data ####
vdem <- readRDS("data/vdem_merged.rds") 

#### Generate Sample ####

vdem <- vdem %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem <- vdem %>%
  mutate(public_goods = v2dlencmps)

vdem_rel <- vdem %>%
  drop_na(public_goods)

vdem_rel <- vdem_rel %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE)

vdem_rel <- vdem_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_rel2 <- vdem_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol, urban_percent)

vdem_rel <- vdem_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, public_goods, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_rel2 <- vdem_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, public_goods, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

### Computing group_mean and de-meaned variables ####

vdem_rel <- sjmisc::de_mean(vdem_rel, pol_incl, v2xps_party, v2x_regime,
                            v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                            grp = country_id)

vdem_rel2 <- sjmisc::de_mean(vdem_rel2, pol_incl, v2xps_party, v2x_regime,
                             e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

#### Descriptive statistics, Appendix E ##

cols1 <- c("public_goods", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols2 <- c("public_goods", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 


# Model 1 and 2 #
stargazer(
  title="Summary Statistics for Public Goods Provision Dataset Dataset", 
  as.data.frame(vdem_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Public Goods Provision", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 3-6 #
stargazer(
  title="Summary Statistics for Public Goods Provision Dataset Dataset", 
  as.data.frame(vdem_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Public Goods Provision", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

## reference Level for Factor ##
vdem_rel$auto_regime_type  <- as.factor(vdem_rel$auto_regime_type)
vdem_rel2$auto_regime_type  <- as.factor(vdem_rel2$auto_regime_type)

vdem_rel$auto_regime_type = relevel(vdem_rel$auto_regime_type, ref=1)
vdem_rel2$auto_regime_type = relevel(vdem_rel2$auto_regime_type, ref=1)


#### Models ####

m1 <- lmer(public_goods ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_rel2, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE

m1.1 <- lmer(public_goods ~  year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               (1  | country_id),
             data = vdem_rel, 
             REML = TRUE, 
             control=lmerControl(optCtrl=list(maxfun=1e9)))
isSingular(m1.1, tol = 1e-05) # FALSE

m2 <- lmer(public_goods ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m2, tol = 1e-05) # FALSE

m2.1 <- lmer(public_goods ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
             data = vdem_rel2)
isSingular(m2.1, tol = 1e-05) # FALSE

# Interaction Models, compare Giesselmann and Schmidt Catran 2018 (double-demeanining of interactions)

m3 <- lmer(public_goods ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m3, tol = 1e-05) # FALSE

m3.1 <- lmer(public_goods ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
               v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
               v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
               (1  | country_id),
             data = vdem_rel2)
isSingular(m3.1, tol = 1e-05) # FALSE

tab_model(m1, m2, m3, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3"))

#### Regression Output ####

texreg(list(m1, m2, m3), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Public Good Spending",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:E3", 
       longtable = TRUE)


texreg(list(m1.1, m2.1, m3.1), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party organizations", 
                             "Party organizations (b)", 
                             "Party Branches", 
                             "Party Branches (b)", 
                             "Party Linkages", 
                             "Party Linkages (b)", 
                             "Distinct Party Platforms", 
                             "Distinct Party Platforms (b)", 
                             "Legislative Party Cohesion", 
                             "Legislative Party Cohesion (b)", 
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * Party organizations", 
                             "CM Autocracy * Party organizations", 
                             "HM Autocracy * Party Branches", 
                             "CM Autocracy * Party Branches", 
                             "HM Autocracy * Party Linkages", 
                             "CM Autocracy * Party Linkages", 
                             "HM Autocracy * Distinct Party Platforms", 
                             "CM Autocracy * Distinct Party Platforms", 
                             "HM Autocracy * Legislative Party Cohesion",
                             "CM Autocracy * Legislative Party Cohesion",
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "Inclusiveness * Party organizations",
                             "Inclusiveness * Party Branches",
                             "Inclusiveness * Party Linkages",
                             "Inclusiveness * Distinct Party Platforms",
                             "Inclusiveness * Legislative Party Cohesion",
                             "Inclusiveness * Party organizations * HM Autocracy",
                             "Inclusiveness * Party organizations * CM Autocracy",
                             "Inclusiveness * Party Branches * HM Autocracy",
                             "Inclusiveness * Party Branches * CM Autocracy",
                             "Inclusiveness * Party Linkages * HM Autocracy",
                             "Inclusiveness * Party Linkages * CM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * HM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * CM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * HM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * CM Autocracy"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Public Good Spending",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:E4", 
       longtable = TRUE)


#### Marginal Plots ####

margin_data <- ggpredict(m2, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "prredicted Public Good Spending") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$public_goods), linetype="dashed")


margin_data <- ggpredict(m2, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted Public Good Spending") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$public_goods), linetype="dashed")


ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 1, nrow = 2, 
          common.legend = TRUE)
ggsave("Output/Margins/S6_Figure_1.pdf", height = 20, width = 24, units= c("cm"))


#### Figure 4 Interaction effects ####

sd(vdem_rel2$pol_incl_dm)

margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "pol_incl_dm[-0.5, 0.5]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.5", "0.5"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted Public Good Spedning") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$public_goods), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)

ggarrange(plot1, 
          ncol = 1, 
          common.legend = TRUE)
ggsave("Output/Margins/S6_Figure_2.pdf", height = 14, width = 24, units= c("cm"))

##########################################################################
##########################################################################

#### APPENDIX F ####
#### 7. OLS Models relative redistribution  ####

# clear workspace
rm(list=ls())

#### Load Data ####
vdem <- readRDS("data/vdem_merged.rds") 

vdem <- vdem %>%
  drop_na(v2x_regime)

#### Generate Sample ####

vdem <- vdem %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_rel <- vdem %>%
  drop_na(rel_red)

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE) # Distinct observations

vdem_rel <- vdem_rel %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_rel <- vdem_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_rel2 <- vdem_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1 )

vdem_rel <- vdem_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, rel_red, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_rel2 <- vdem_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, rel_red, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

## reference Level for Factor ##
vdem_rel$auto_regime_type  <- as.factor(vdem_rel$auto_regime_type)
vdem_rel2$auto_regime_type  <- as.factor(vdem_rel2$auto_regime_type)

vdem_rel$auto_regime_type = relevel(vdem_rel$auto_regime_type, ref=1)
vdem_rel2$auto_regime_type = relevel(vdem_rel2$auto_regime_type, ref=1)

#### Model with one year time lag and fixed (country and time) structure ####

library(plm)

mod.1 = plm(rel_red ~  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel, 
            model="within", 
            effect = "twoways")
summary(mod.1)

mod.2 = plm(rel_red ~  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="within", 
            effect = "twoways")
summary(mod.2)

mod.3 = plm(rel_red ~  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(v2xps_party, 1)* plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + 
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="within", 
            effect = "twoways")
summary(mod.3)

# Main models

mods = list(mod.1, mod.2, mod.3)
vcm = lapply(mods, function(x) plm::vcovBK(x, cluster="time", pairwise = FALSE)) # # Beck-Katz SEs for unbalanced panel
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mods)) {pval[[i]] = pt(abs(coef(mods[[i]]) / rse[[i]]), df=1543, lower.tail=FALSE) * 2}
sapply(mods, function(x) round(sigma(x), 3))
sapply(mods, function(x) round(pbgtest(x, order=1)$p.value, 3)) # Wooldridge test for serial correlation

texreg(mods,
       override.se=rse, override.pvalues=pval, 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("Pol Inclusiveness",
                             "Party Institutionalization",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "Pol Inclusiveness * HM Autocracy", 
                             "Pol Inclusiveness * CM Autocracy", 
                             "Party Institutionalization * HM Autocracy", 
                             "Party Institutionalization * CM Autocracy", 
                             "GDP pc log",
                             "Population log",
                             "Communist Ideology", 
                             "Urban Percentage", 
                             "Manufacturing Sector", 
                             "Party Institutionalization * Inclusiveness",
                             "HM Autocracy *Party Institutionalization * Inclusiveness", 
                             "CM Autocracy * Party Institutionalization * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Fixed Effects Model prediciting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:F1", 
       longtable = TRUE)

##### Table F2: including Lagged DV ####

mod.1 = plm(rel_red ~  plm::lag(rel_red, 1) + plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel, 
            model="within", 
            effect = "twoways")
summary(mod.1)

mod.2 = plm(rel_red ~ plm::lag(rel_red, 1) + plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="within", 
            effect = "twoways")
summary(mod.2)

mod.3 = plm(rel_red ~  plm::lag(rel_red, 1)  + plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(v2xps_party, 1)* plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + 
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="within", 
            effect = "twoways")
summary(mod.3)


# Main models

mods = list(mod.1, mod.2, mod.3)
vcm = lapply(mods, function(x) plm::vcovBK(x, cluster="time", pairwise = FALSE)) # # Beck-Katz SEs for unbalanced panel
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mods)) {pval[[i]] = pt(abs(coef(mods[[i]]) / rse[[i]]), df=1542, lower.tail=FALSE) * 2}
sapply(mods, function(x) round(sigma(x), 3))
sapply(mods, function(x) round(pbgtest(x, order=1)$p.value, 3)) # Wooldridge test for serial correlation

texreg(mods,
       override.se=rse, override.pvalues=pval, 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("Lagged DV", 
                             "Pol Inclusiveness",
                             "Party Institutionalization",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "Pol Inclusiveness * HM Autocracy", 
                             "Pol Inclusiveness * CM Autocracy", 
                             "Party Institutionalization * HM Autocracy", 
                             "Party Institutionalization * CM Autocracy", 
                             "GDP pc log",
                             "Population log",
                             "Communist Ideology", 
                             "Urban Percentage", 
                             "Manufacturing Sector", 
                             "Party Institutionalization * Inclusiveness",
                             "HM Autocracy *Party Institutionalization * Inclusiveness", 
                             "CM Autocracy * Party Institutionalization * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Fixed Effects Model prediciting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:F2", 
       longtable = TRUE)


#### Model with one year time lag and Pooled stracture ####
## F3 ##


mod.1 = plm(rel_red ~  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel, 
            model="pooling", 
            effect = "twoways")
summary(mod.1)

mod.2 = plm(rel_red ~  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="pooling", 
            effect = "twoways")
summary(mod.2)

mod.3 = plm(rel_red ~  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(v2xps_party, 1)* plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + 
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="pooling", 
            effect = "twoways")
summary(mod.3)


# Main models

mods = list(mod.1, mod.2, mod.3)
vcm = lapply(mods, function(x) plm::vcovBK(x, cluster="time", pairwise = FALSE)) # # Beck-Katz SEs for unbalanced panel
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mods)) {pval[[i]] = pt(abs(coef(mods[[i]]) / rse[[i]]), df=1692, lower.tail=FALSE) * 2}
sapply(mods, function(x) round(sigma(x), 3))
sapply(mods, function(x) round(pbgtest(x, order=1)$p.value, 3)) # Wooldridge test for serial correlation

texreg(mods,
       override.se=rse, override.pvalues=pval, 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("Intercept",
                             "Pol Inclusiveness",
                             "Party Institutionalization",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "Pol Inclusiveness * HM Autocracy", 
                             "Pol Inclusiveness * CM Autocracy", 
                             "Party Institutionalization * HM Autocracy", 
                             "Party Institutionalization * CM Autocracy", 
                             "GDP pc log",
                             "Population log",
                             "Communist Ideology", 
                             "Urban Percentage", 
                             "Manufacturing Sector", 
                             "Party Institutionalization * Inclusiveness",
                             "HM Autocracy *Party Institutionalization * Inclusiveness", 
                             "CM Autocracy * Party Institutionalization * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Fixed Effects Model prediciting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:F3", 
       longtable = TRUE)

#### Model with one year time lag and DV as Control and Pooled structure #####


mod.1 = plm(rel_red ~  plm::lag(rel_red, 1) + plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel, 
            model="pooling", 
            effect = "twoways")
summary(mod.1)

mod.2 = plm(rel_red ~  plm::lag(rel_red, 1) + plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="pooling", 
            effect = "twoways")
summary(mod.2)

mod.3 = plm(rel_red ~  plm::lag(rel_red, 1) +  plm::lag(pol_incl, 1) + plm::lag(v2xps_party, 1) + plm::lag(auto_regime_type, 1) + 
              plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + plm::lag(v2xps_party, 1)*plm::lag(auto_regime_type, 1) +
              plm::lag(v2xps_party, 1)* plm::lag(pol_incl, 1)* plm::lag(auto_regime_type, 1) + 
              plm::lag(e_migdppcln, 1) + plm::lag(e_wb_pop_ln, 1) + 
              plm::lag(v2exl_legitideolcr_1, 1) + plm::lag(urban_percent, 1) + plm::lag(manufacturing_percent_ipol, 1), 
            index=c("country_id", "year"), 
            data = vdem_rel2, 
            model="pooling", 
            effect = "twoways")
summary(mod.3)


# Main models

mods = list(mod.1, mod.2, mod.3)
vcm = lapply(mods, function(x) plm::vcovBK(x, cluster="time", pairwise = FALSE)) # # Beck-Katz SEs for unbalanced panel
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mods)) {pval[[i]] = pt(abs(coef(mods[[i]]) / rse[[i]]), df=1691, lower.tail=FALSE) * 2}
sapply(mods, function(x) round(sigma(x), 3))
sapply(mods, function(x) round(pbgtest(x, order=1)$p.value, 3)) # Wooldridge test for serial correlation

texreg(mods,
       override.se=rse, override.pvalues=pval, 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("Intercept",
                             "Lagged DV", 
                             "Pol Inclusiveness",
                             "Party Institutionalization",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "Pol Inclusiveness * HM Autocracy", 
                             "Pol Inclusiveness * CM Autocracy", 
                             "Party Institutionalization * HM Autocracy", 
                             "Party Institutionalization * CM Autocracy", 
                             "GDP pc log",
                             "Population log",
                             "Communist Ideology", 
                             "Urban Percentage", 
                             "Manufacturing Sector", 
                             "Party Institutionalization * Inclusiveness",
                             "HM Autocracy *Party Institutionalization * Inclusiveness", 
                             "CM Autocracy * Party Institutionalization * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Fixed Effects Model prediciting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:F4", 
       longtable = TRUE)

